function [yExtra_cu,labels,yExtraSS] = ProjectExtraControls(y_cu,x_cu,model)

% The required parameters
params  = model.params;
DELTA   = params.DELTA;
ZI      = params.ZI;
KAPAw   = params.KAPAw;
PAIss   = params.PAIss;
NU      = params.NU;
Kss     = params.Kss;
lss     = params.lss;
THETA   = params.THETA;
PHIzero = params.PHIzero;
PHI     = params.PHI;
BETTA   = params.BETTA;
ALFA    = params.ALFA;
B       = params.B;
CHI0    = params.CHI0;
CHI     = params.CHI;
MUZss   = params.MUZss;
U0      = params.U0;
U0d     = params.U0d;
Css     = model.params.Css;
OUTPUTss= model.params.OUTPUTss;
Wss     = model.params.Wss;
Rss     = model.params.Rss;
AA      = model.params.AA;

% The required variables on the considered grid
pos_c   = find(strcmp(model.labely,'$c_t$'));
pos_pai = find(strcmp(model.labely,'$\pi_t$'));
pos_evf = find(strcmp(model.labely,'$evf_t$'));
c_cu    = exp(y_cu(pos_c,:));
pai_cu  = exp(y_cu(pos_pai,:));
evf_cu  = exp(y_cu(pos_evf,:));

c_ba1      = exp(x_cu(1,:)+model.hSS(1,1));
muz_cu     = exp(x_cu(2,:)+model.hSS(2,1));
d_cu       = exp(x_cu(3,:)+model.hSS(3,1));
n_cu       = exp(x_cu(4,:)+model.hSS(4,1));
paistar_cu = exp(x_cu(5,:)+model.hSS(5,1));
a_cu       = exp(x_cu(6,:)+model.hSS(6,1));


%EQ7: Agg. ressource constraint
output_cu  = (c_cu  + DELTA*Kss)./(1-ZI/2.*(pai_cu ./PAIss.^NU-1).^2);

% EQ 8: The aggregated production function to get output. 
l_cu  = (output_cu ./(a_cu *Kss^THETA)).^(1/(1-THETA));

% EQ 2: Household's First-order condition for labour
w_cu = KAPAw*Wss+(1-KAPAw)*(n_cu.*PHIzero.*(1-l_cu).^(-1/PHI)./((c_cu -B*c_ba1./muz_cu)/Css^CHI0).^(-CHI)*Css^CHI0);

% EQ 4: Firms FOC for labour
mc_cu  = w_cu./((1-THETA).*a_cu.*Kss^THETA.*l_cu.^-THETA);
MCss   = Wss /((1-THETA).*1.*Kss^THETA.*lss.^-THETA);

% Dividends
div_cu = c_cu - w_cu*l_cu;

% The value function
vf_cu = - (U0 + d_cu.*(1/(1-CHI).*((c_cu -B*c_ba1./muz_cu)./Css^CHI0).^(1-CHI)+U0d ...
                      + n_cu*PHIzero./(1-1/PHI).*(1-l_cu ).^(1-1/PHI)) ) ...
        + BETTA.*AA.*evf_cu.^(1/(1-ALFA));
VFss  = -(U0+U0d+1/(1-CHI)*((Css-B*Css*MUZss^-1)/Css^CHI0)^(1-CHI) + PHIzero*(1-lss)^(1-1/PHI)/(1-1/PHI))/(1-BETTA*MUZss^((1-CHI)*(1-CHI0)));
%% Saving 
% Getting the additional controls log transformed
yExtra_cu.output_cu = log(output_cu);
yExtra_cu.l_cu      = log(l_cu);
yExtra_cu.w_cu      = log(w_cu);
yExtra_cu.mc_cu     = log(mc_cu);
yExtra_cu.div_cu    = log(max(div_cu,1D-4));  
yExtra_cu.vf_cu     = log(vf_cu);

labels = {'$output_t$','$l_t$','$w_t$','$mc_t$','$div_t$','$vf_t$'};
yExtraSS.output_cu = log(OUTPUTss);
yExtraSS.l_cu      = log(model.params.lss);
yExtraSS.w_cu      = log(model.params.Wss);
yExtraSS.mc_cu     = log(MCss);
yExtraSS.div_cu    = log(Css-model.params.Wss*model.params.lss);
yExtraSS.vf_cu     = log(VFss);
end
